setwd("[your home directory]/IRF_Figures/InputData")
rm(list = ls())
library(tidyselect)
library(readr)
library(lpirfs)
library(gridExtra)
library(ggpubr)
library(cowplot)
library(colorspace)
library(foreign)
library(haven)
library(magrittr)
library(dplyr)

#-----------------------------Clean Data For Charts-------------------------####
mindate=1
maxdate=121 # 121= run through 2019, 92=run thru 2012q2 


dfa_gini=read_csv("dfa_gini_networth_low.csv")
dfa_networth_shares <- read_csv("dfa_shares_bus_cycle2_networth_2019q3.csv") 

CE_shares=dfa_networth_shares['corp_equ_mfs']
RE_shares=dfa_networth_shares['real_estate']
NW_shares=dfa_networth_shares['networth']
ASSETS_shares=dfa_networth_shares['assets']
Group_shares=t(matrix(unlist(NW_shares), nrow=6))
Group_shares=Group_shares[mindate:maxdate,]
Group_shares2=data.frame(rowSums(Group_shares[,1:2]), Group_shares[,3], rowSums(Group_shares[,4:5]))

Group_shares3=data.frame(rowSums(Group_shares[,1:3]),   rowSums(Group_shares[,4:6]))

TOP01=data.frame(Group_shares[,1:1])
TOP1a=data.frame(rowSums(Group_shares[,1:2]))
TOP10a=data.frame(rowSums(Group_shares[,1:3]))
TOP1=data.frame(Group_shares[,2:2])
TOP10=data.frame(Group_shares[,3:3])
BOT90=data.frame(rowSums(Group_shares[,4:6]))
P7090=data.frame(Group_shares[,4:4])
P5070=data.frame(Group_shares[,5:5])
BOT50=data.frame(Group_shares[,6:6])

Group_shares4=data.frame(TOP01,TOP1,TOP10)
colnames(Group_shares4) <- c("TOP01","TOP1","TOP10")
Group_shares5=data.frame(TOP01,TOP1)
colnames(Group_shares5) <- c("TOP01","TOP1")

CE_Group_shares=t(matrix(unlist(CE_shares), nrow=6))
CE_Group_shares=CE_Group_shares[mindate:maxdate,]
CE_Group_shares2=data.frame(rowSums(CE_Group_shares[,1:2]), CE_Group_shares[,3], rowSums(CE_Group_shares[,4:5]))


RE_Group_shares=t(matrix(unlist(RE_shares), nrow=6))
RE_Group_shares=RE_Group_shares[mindate:maxdate,]
RE_Group_shares2=data.frame(rowSums(RE_Group_shares[,1:2]), RE_Group_shares[,3], rowSums(RE_Group_shares[,4:5]))


ASSETS_Group_shares=t(matrix(unlist(ASSETS_shares), nrow=6))
ASSETS_Group_shares=ASSETS_Group_shares[mindate:maxdate,]
ASSETS_Group_shares2=data.frame(rowSums(ASSETS_Group_shares[,1:2]), ASSETS_Group_shares[,3], rowSums(ASSETS_Group_shares[,4:6]))



plotnames=c('toppoint1','top1', 'next9','prctile60to90','prctile30to60','prctile0to30')
plotnames2=c('top1', 'next9','prctile10to100')
plotnames3=c('top10','prctile10to100')

DFAGINI=dfa_gini['gini']
GINI=DFAGINI[mindate:maxdate,]

FFrate <- read.csv("FEDFUNDS.csv")
FFrate=FFrate[mindate:maxdate,]
FFR=FFrate['FEDFUNDS']
FFR2=c(9.73,FFR[,1])
FFR_diff=diff(FFR2,1)
FFR_diff=c(FFR_diff)


housing <- read.csv("CSUSHPISA_June2023.csv")
#annualized ror
house=housing['annualized']
HOUSE=as.data.frame(house[mindate:maxdate,])

# upload_shocks_june2023.csv -- includes the Gertler and Karadi 2015 shocks
MPShocks_all<- read.csv("upload_shocks_june2023.csv")
MPS_RR_shocks=MPShocks_all['RR_Original']
MPS_RR_shocks_ext=MPShocks_all['RR_Extended']
MPS_NS_shocks1=MPShocks_all['NS_FFR']
MPS_NS_shocks2=MPShocks_all['NS_News']
MPS_BRW_shocks=MPShocks_all['BRW']


# Romer and Romer
MPS_RR=as.data.frame(MPS_RR_shocks[mindate:maxdate,])
MPS_RR_ext=as.data.frame(MPS_RR_shocks_ext[mindate:maxdate,])
# Nakamura Steinsson
MPS_NS1=as.data.frame(MPS_NS_shocks1[mindate:maxdate,])
MPS_NS2=as.data.frame(MPS_NS_shocks2[mindate:maxdate,])
MPS_BRW=as.data.frame(MPS_BRW_shocks[mindate:maxdate,])

# New Gertler shocks (June 2023)
MPS_GK_shock1=MPShocks_all['GK_MP1']
MPS_GK_shock2=MPShocks_all['GK_FF4']

MPS_GK1=as.data.frame(MPS_GK_shock1[mindate:maxdate,])
MPS_GK2=as.data.frame(MPS_GK_shock2[mindate:maxdate,])


Inflation3<- read.csv("CPI.csv")
Inflation3=Inflation3[mindate:maxdate,]
CPI=Inflation3['CPIAUCSL_PCH']
INF=CPI

gdp<- read.csv("GDPGAP.csv")
gdp=gdp[mindate:maxdate,]
GDP=gdp['GDPC1_GDPPOT']

unemployment<- read.csv("UNRATE.csv")
unemployment=unemployment[mindate:maxdate,]
UNEMP=unemployment['UNRATE']

ugap<- read.csv("ugap.csv")
ugap=ugap[mindate:maxdate,]
UGAP=ugap['NROU_UNRATE']

sp500 <- read.csv("upload_sp500.csv")
sp500 =sp500['Return_ann']

sp=sp500[(mindate:maxdate),]
SP500=sp
STOCK= sp


plotselect=seq(6,36,6)
plotselect2=seq(4,16,4)
plotselect3=seq(3,9,3)

#--------------------HOUSE STOCK INF UNEMP GDP MP Top .1%-------------------####
shock="ALL7x7_top01"

IRF_data_top016=as.data.frame(cbind(TOP01,HOUSE,STOCK,INF,UNEMP,GDP,-FFR))
results_lin_top016= lp_lin(IRF_data_top016, 
                           lags_endog_lin = 8,    # Number of lags for endogenous data
                           trend          = 1,    # 0 = no trend, 1 = trend, 2 = trend & trend^2    
                           shock_type     = 0,    # 0 = standard deviation shock, 1 = unit shock
                           confint        = 1.645,
                           hor=16
)

# save all plots for top 0.1% wealth share IRF
linear_plots_top01_all7 <- plot_lin(results_lin_top016)
# save each plot individually (grab top row of the 7x7, one-by-one)
p1_01 <- linear_plots_top01_all7[[1]]
p2_01 <- linear_plots_top01_all7[[2]]
p3_01 <- linear_plots_top01_all7[[3]]
p4_01 <- linear_plots_top01_all7[[4]]
p5_01 <- linear_plots_top01_all7[[5]]
p6_01 <- linear_plots_top01_all7[[6]]
p7_01 <- linear_plots_top01_all7[[7]]

#--------------------HOUSE STOCK INF UNEMP GDP MP Top 99-99.9%--------------####
shock="ALL7x7_top1"

IRF_data_top16=as.data.frame(cbind(TOP1,HOUSE,STOCK,INF,UNEMP,GDP,-FFR))
results_lin_top16= lp_lin(IRF_data_top16, 
                          lags_endog_lin = 8,    # Number of lags for endogenous data
                          trend          = 1,    # 0 = no trend, 1 = trend, 2 = trend & trend^2    
                          shock_type     = 0,    # 0 = standard deviation shock, 1 = unit shock
                          confint        = 1.645,
                          hor=16
)
linear_plots_top1_all7 <- plot_lin(results_lin_top16)

p1_1 <- linear_plots_top1_all7[[1]]
p2_1 <- linear_plots_top1_all7[[2]]
p3_1 <- linear_plots_top1_all7[[3]]
p4_1 <- linear_plots_top1_all7[[4]]
p5_1 <- linear_plots_top1_all7[[5]]
p6_1 <- linear_plots_top1_all7[[6]]
p7_1 <- linear_plots_top1_all7[[7]]

#----------------HOUSE STOCK INF UNEMP GDP MP Top 90-99%--------------------####
shock="ALL7x7_top10"

IRF_data_top106=as.data.frame(cbind(TOP10,HOUSE,STOCK,INF,UNEMP,GDP,-FFR))
results_lin_top106= lp_lin(IRF_data_top106, 
                           lags_endog_lin = 8,    # Number of lags for endogenous data
                           trend          = 1,    # 0 = no trend, 1 = trend, 2 = trend & trend^2    
                           shock_type     = 0,    # 0 = standard deviation shock, 1 = unit shock
                           confint        = 1.645,
                           hor=16
)
linear_plots_top10_all7 <- plot_lin(results_lin_top106)

p1_10 <- linear_plots_top10_all7[[1]]
p2_10 <- linear_plots_top10_all7[[2]]
p3_10 <- linear_plots_top10_all7[[3]]
p4_10 <- linear_plots_top10_all7[[4]]
p5_10 <- linear_plots_top10_all7[[5]]
p6_10 <- linear_plots_top10_all7[[6]]
p7_10 <- linear_plots_top10_all7[[7]]

#-------------------HOUSE STOCK INF UNEMP GDP MP Top 70-90%-----------------####
shock="ALL7x7_p7090"

IRF_data_p70906=as.data.frame(cbind(P7090,HOUSE,STOCK,INF,UNEMP,GDP,-FFR))
results_lin_p70906= lp_lin(IRF_data_p70906, 
                           lags_endog_lin = 8,    # Number of lags for endogenous data
                           trend          = 1,    # 0 = no trend, 1 = trend, 2 = trend & trend^2    
                           shock_type     = 0,    # 0 = standard deviation shock, 1 = unit shock
                           confint        = 1.645,
                           hor=16
)
linear_plots_p7090_all7 <- plot_lin(results_lin_p70906)

p1_79 <- linear_plots_p7090_all7[[1]]
p2_79 <- linear_plots_p7090_all7[[2]]
p3_79 <- linear_plots_p7090_all7[[3]]
p4_79 <- linear_plots_p7090_all7[[4]]
p5_79 <- linear_plots_p7090_all7[[5]]
p6_79 <- linear_plots_p7090_all7[[6]]
p7_79 <- linear_plots_p7090_all7[[7]]

#----------------HOUSE STOCK INF UNEMP GDP MP Top 50-70%--------------------####
shock="ALL7x7_p5070"

IRF_data_p50706=as.data.frame(cbind(P5070,HOUSE,STOCK,INF,UNEMP,GDP,-FFR))
results_lin_p50706= lp_lin(IRF_data_p50706, 
                           lags_endog_lin = 8,    # Number of lags for endogenous data
                           trend          = 1,    # 0 = no trend, 1 = trend, 2 = trend & trend^2    
                           shock_type     = 0,    # 0 = standard deviation shock, 1 = unit shock
                           confint        = 1.645,
                           hor=16
)
linear_plots_p5070_all7 <- plot_lin(results_lin_p50706)

p1_57 <- linear_plots_p5070_all7[[1]]
p2_57 <- linear_plots_p5070_all7[[2]]
p3_57 <- linear_plots_p5070_all7[[3]]
p4_57 <- linear_plots_p5070_all7[[4]]
p5_57 <- linear_plots_p5070_all7[[5]]
p6_57 <- linear_plots_p5070_all7[[6]]
p7_57 <- linear_plots_p5070_all7[[7]]

#------------------HOUSE STOCK INF UNEMP GDP MP Bottom 50%------------------####
shock="ALL7x7_bot50"

IRF_data_bot506=as.data.frame(cbind(BOT50,HOUSE,STOCK,INF,UNEMP,GDP,-FFR))
results_lin_bot506= lp_lin(IRF_data_bot506, 
                           lags_endog_lin = 8,    # Number of lags for endogenous data
                           trend          = 1,    # 0 = no trend, 1 = trend, 2 = trend & trend^2    
                           shock_type     = 0,    # 0 = standard deviation shock, 1 = unit shock
                           confint        = 1.645,
                           hor=16
)
linear_plots_bot50_all7 <- plot_lin(results_lin_bot506)

p1_b5 <- linear_plots_bot50_all7[[1]]
p2_b5 <- linear_plots_bot50_all7[[2]]
p3_b5 <- linear_plots_bot50_all7[[3]]
p4_b5 <- linear_plots_bot50_all7[[4]]
p5_b5 <- linear_plots_bot50_all7[[5]]
p6_b5 <- linear_plots_bot50_all7[[6]]
p7_b5 <- linear_plots_bot50_all7[[7]]
#----------------------------------Export Series----------------------------####
setwd("[your home directory]/IRF_Figures/OutputData/")
export_figure18_series <- function()
{
  # Create Figure 18a Series 
  figure18aseries_high <- as.data.frame(t(results_lin_top016$irf_lin_up[, , 2]))
  figure18aseries_high <- as.data.frame(figure18aseries_high$V1)
  
  figure18aseries_mean <- as.data.frame(t(results_lin_top016$irf_lin_mean[, , 2]))
  figure18aseries_mean <- as.data.frame(figure18aseries_mean$V1)
  
  figure18aseries_low <- as.data.frame(t(results_lin_top016$irf_lin_low[, , 2]))
  figure18aseries_low <- as.data.frame(figure18aseries_low$V1)
  
  # Create Figure 18b Series 
  figure18bseries_high <- as.data.frame(t(results_lin_top16$irf_lin_up[, , 2]))
  figure18bseries_high <- as.data.frame(figure18bseries_high$V1)
  
  figure18bseries_mean <- as.data.frame(t(results_lin_top16$irf_lin_mean[, , 2]))
  figure18bseries_mean <- as.data.frame(figure18bseries_mean$V1)
  
  figure18bseries_low <- as.data.frame(t(results_lin_top16$irf_lin_low[, , 2]))
  figure18bseries_low <- as.data.frame(figure18bseries_low$V1)
  
  # Create Figure 18c Series 
  figure18cseries_high <- as.data.frame(t(results_lin_top106$irf_lin_up[, , 2]))
  figure18cseries_high <- as.data.frame(figure18cseries_high$V1)
  
  figure18cseries_mean <- as.data.frame(t(results_lin_top106$irf_lin_mean[, , 2]))
  figure18cseries_mean <- as.data.frame(figure18cseries_mean$V1)
  
  figure18cseries_low <- as.data.frame(t(results_lin_top106$irf_lin_low[, , 2]))
  figure18cseries_low <- as.data.frame(figure18cseries_low$V1)
  
  # Create Figure 18d Series 
  figure18dseries_high <- as.data.frame(t(results_lin_p70906$irf_lin_up[, , 2]))
  figure18dseries_high <- as.data.frame(figure18dseries_high$V1)
  
  figure18dseries_mean <- as.data.frame(t(results_lin_p70906$irf_lin_mean[, , 2]))
  figure18dseries_mean <- as.data.frame(figure18dseries_mean$V1)
  
  figure18dseries_low <- as.data.frame(t(results_lin_p70906$irf_lin_low[, , 2]))
  figure18dseries_low <- as.data.frame(figure18dseries_low$V1)
  
  # Create Figure 18e Series 
  figure18eseries_high <- as.data.frame(t(results_lin_p50706$irf_lin_up[, , 2]))
  figure18eseries_high <- as.data.frame(figure18eseries_high$V1)
  
  figure18eseries_mean <- as.data.frame(t(results_lin_p50706$irf_lin_mean[, , 2]))
  figure18eseries_mean <- as.data.frame(figure18eseries_mean$V1)
  
  figure18eseries_low <- as.data.frame(t(results_lin_p50706$irf_lin_low[, , 2]))
  figure18eseries_low <- as.data.frame(figure18eseries_low$V1)
  
  # Create Figure 18f Series 
  figure18fseries_high <- as.data.frame(t(results_lin_bot506$irf_lin_up[, , 2]))
  figure18fseries_high <- as.data.frame(figure18fseries_high$V1)
  
  figure18fseries_mean <- as.data.frame(t(results_lin_bot506$irf_lin_mean[, , 2]))
  figure18fseries_mean <- as.data.frame(figure18fseries_mean$V1)
  
  figure18fseries_low <- as.data.frame(t(results_lin_bot506$irf_lin_low[, , 2]))
  figure18fseries_low <- as.data.frame(figure18fseries_low$V1)
  
  # Merge all series together
  figure18series <- list(figure18aseries_high, figure18aseries_mean, figure18aseries_low, figure18bseries_high, figure18bseries_mean, figure18bseries_low, figure18cseries_high, figure18cseries_mean, figure18cseries_low, figure18dseries_high, figure18dseries_mean, figure18dseries_low, figure18eseries_high, figure18eseries_mean, figure18eseries_low, figure18fseries_high, figure18fseries_mean, figure18fseries_low)
  figure18series <- as.data.frame(figure18series)
  figure18series <- (mutate(figure18series, quarter=row_number()))
  quarter <- as.data.frame(figure18series$quarter)
  figure18series <- list(quarter, figure18aseries_high, figure18aseries_mean, figure18aseries_low, figure18bseries_high, figure18bseries_mean, figure18bseries_low, figure18cseries_high, figure18cseries_mean, figure18cseries_low, figure18dseries_high, figure18dseries_mean, figure18dseries_low, figure18eseries_high, figure18eseries_mean, figure18eseries_low, figure18fseries_high, figure18fseries_mean, figure18fseries_low)
  figure18series <- as.data.frame(figure18series)
  figure18series <- rename(figure18series, quarters = `figure18series.quarter`, figure18ahigh = `figure18aseries_high.V1`, figure18bhigh = `figure18bseries_high.V1`, figure18chigh = `figure18cseries_high.V1`, figure18dhigh = `figure18dseries_high.V1`, figure18ehigh = `figure18eseries_high.V1`, figure18fhigh = `figure18fseries_high.V1`, figure18amean = `figure18aseries_mean.V1`, figure18bmean = `figure18bseries_mean.V1`, figure18cmean = `figure18cseries_mean.V1`, figure18dmean = `figure18dseries_mean.V1`, figure18emean = `figure18eseries_mean.V1`, figure18fmean = `figure18fseries_mean.V1`, figure18alow = `figure18aseries_low.V1`, figure18blow = `figure18bseries_low.V1`, figure18clow = `figure18cseries_low.V1`, figure18dlow = `figure18dseries_low.V1`, figure18elow = `figure18eseries_low.V1`, figure18flow = `figure18fseries_low.V1`)
  
  # Export the series 
  write.csv(figure18series, "figure18series.csv")
}
export_figure18_series()

export_figure19_series <- function()
{
  # Create Figure 19a Series 
  figure19aseries_high <- as.data.frame(t(results_lin_top016$irf_lin_up[, , 3]))
  figure19aseries_high <- as.data.frame(figure19aseries_high$V1)
  
  figure19aseries_mean <- as.data.frame(t(results_lin_top016$irf_lin_mean[, , 3]))
  figure19aseries_mean <- as.data.frame(figure19aseries_mean$V1)
  
  figure19aseries_low <- as.data.frame(t(results_lin_top016$irf_lin_low[, , 3]))
  figure19aseries_low <- as.data.frame(figure19aseries_low$V1)
  
  # Create Figure 19b Series 
  figure19bseries_high <- as.data.frame(t(results_lin_top16$irf_lin_up[, , 3]))
  figure19bseries_high <- as.data.frame(figure19bseries_high$V1)
  
  figure19bseries_mean <- as.data.frame(t(results_lin_top16$irf_lin_mean[, , 3]))
  figure19bseries_mean <- as.data.frame(figure19bseries_mean$V1)
  
  figure19bseries_low <- as.data.frame(t(results_lin_top16$irf_lin_low[, , 3]))
  figure19bseries_low <- as.data.frame(figure19bseries_low$V1)
  
  # Create Figure 19c Series 
  figure19cseries_high <- as.data.frame(t(results_lin_top106$irf_lin_up[, , 3]))
  figure19cseries_high <- as.data.frame(figure19cseries_high$V1)
  
  figure19cseries_mean <- as.data.frame(t(results_lin_top106$irf_lin_mean[, , 3]))
  figure19cseries_mean <- as.data.frame(figure19cseries_mean$V1)
  
  figure19cseries_low <- as.data.frame(t(results_lin_top106$irf_lin_low[, , 3]))
  figure19cseries_low <- as.data.frame(figure19cseries_low$V1)
  
  # Create Figure 19d Series 
  figure19dseries_high <- as.data.frame(t(results_lin_p70906$irf_lin_up[, , 3]))
  figure19dseries_high <- as.data.frame(figure19dseries_high$V1)
  
  figure19dseries_mean <- as.data.frame(t(results_lin_p70906$irf_lin_mean[, , 3]))
  figure19dseries_mean <- as.data.frame(figure19dseries_mean$V1)
  
  figure19dseries_low <- as.data.frame(t(results_lin_p70906$irf_lin_low[, , 3]))
  figure19dseries_low <- as.data.frame(figure19dseries_low$V1)
  
  # Create Figure 19e Series 
  figure19eseries_high <- as.data.frame(t(results_lin_p50706$irf_lin_up[, , 3]))
  figure19eseries_high <- as.data.frame(figure19eseries_high$V1)
  
  figure19eseries_mean <- as.data.frame(t(results_lin_p50706$irf_lin_mean[, , 3]))
  figure19eseries_mean <- as.data.frame(figure19eseries_mean$V1)
  
  figure19eseries_low <- as.data.frame(t(results_lin_p50706$irf_lin_low[, , 3]))
  figure19eseries_low <- as.data.frame(figure19eseries_low$V1)
  
  # Create Figure 19f Series 
  figure19fseries_high <- as.data.frame(t(results_lin_bot506$irf_lin_up[, , 3]))
  figure19fseries_high <- as.data.frame(figure19fseries_high$V1)
  
  figure19fseries_mean <- as.data.frame(t(results_lin_bot506$irf_lin_mean[, , 3]))
  figure19fseries_mean <- as.data.frame(figure19fseries_mean$V1)
  
  figure19fseries_low <- as.data.frame(t(results_lin_bot506$irf_lin_low[, , 3]))
  figure19fseries_low <- as.data.frame(figure19fseries_low$V1)
  
  # Merge all series together
  figure19series <- list(figure19aseries_high, figure19aseries_mean, figure19aseries_low, figure19bseries_high, figure19bseries_mean, figure19bseries_low, figure19cseries_high, figure19cseries_mean, figure19cseries_low, figure19dseries_high, figure19dseries_mean, figure19dseries_low, figure19eseries_high, figure19eseries_mean, figure19eseries_low, figure19fseries_high, figure19fseries_mean, figure19fseries_low)
  figure19series <- as.data.frame(figure19series)
  figure19series <- (mutate(figure19series, quarter=row_number()))
  quarter <- as.data.frame(figure19series$quarter)
  figure19series <- list(quarter, figure19aseries_high, figure19aseries_mean, figure19aseries_low, figure19bseries_high, figure19bseries_mean, figure19bseries_low, figure19cseries_high, figure19cseries_mean, figure19cseries_low, figure19dseries_high, figure19dseries_mean, figure19dseries_low, figure19eseries_high, figure19eseries_mean, figure19eseries_low, figure19fseries_high, figure19fseries_mean, figure19fseries_low)
  figure19series <- as.data.frame(figure19series)
  figure19series <- rename(figure19series, quarters = `figure19series.quarter`, figure19ahigh = `figure19aseries_high.V1`, figure19bhigh = `figure19bseries_high.V1`, figure19chigh = `figure19cseries_high.V1`, figure19dhigh = `figure19dseries_high.V1`, figure19ehigh = `figure19eseries_high.V1`, figure19fhigh = `figure19fseries_high.V1`, figure19amean = `figure19aseries_mean.V1`, figure19bmean = `figure19bseries_mean.V1`, figure19cmean = `figure19cseries_mean.V1`, figure19dmean = `figure19dseries_mean.V1`, figure19emean = `figure19eseries_mean.V1`, figure19fmean = `figure19fseries_mean.V1`, figure19alow = `figure19aseries_low.V1`, figure19blow = `figure19bseries_low.V1`, figure19clow = `figure19cseries_low.V1`, figure19dlow = `figure19dseries_low.V1`, figure19elow = `figure19eseries_low.V1`, figure19flow = `figure19fseries_low.V1`)
  
  # Export the series 
  write.csv(figure19series, "figure19series.csv")
}
export_figure19_series()

export_figure20_series <- function()
{
  # Create Figure 20a Series 
  figure20aseries_high <- as.data.frame(t(results_lin_top016$irf_lin_up[, , 4]))
  figure20aseries_high <- as.data.frame(figure20aseries_high$V1)
  
  figure20aseries_mean <- as.data.frame(t(results_lin_top016$irf_lin_mean[, , 4]))
  figure20aseries_mean <- as.data.frame(figure20aseries_mean$V1)
  
  figure20aseries_low <- as.data.frame(t(results_lin_top016$irf_lin_low[, , 4]))
  figure20aseries_low <- as.data.frame(figure20aseries_low$V1)
  
  # Create Figure 20b Series 
  figure20bseries_high <- as.data.frame(t(results_lin_top16$irf_lin_up[, , 4]))
  figure20bseries_high <- as.data.frame(figure20bseries_high$V1)
  
  figure20bseries_mean <- as.data.frame(t(results_lin_top16$irf_lin_mean[, , 4]))
  figure20bseries_mean <- as.data.frame(figure20bseries_mean$V1)
  
  figure20bseries_low <- as.data.frame(t(results_lin_top16$irf_lin_low[, , 4]))
  figure20bseries_low <- as.data.frame(figure20bseries_low$V1)
  
  # Create Figure 20c Series 
  figure20cseries_high <- as.data.frame(t(results_lin_top106$irf_lin_up[, , 4]))
  figure20cseries_high <- as.data.frame(figure20cseries_high$V1)
  
  figure20cseries_mean <- as.data.frame(t(results_lin_top106$irf_lin_mean[, , 4]))
  figure20cseries_mean <- as.data.frame(figure20cseries_mean$V1)
  
  figure20cseries_low <- as.data.frame(t(results_lin_top106$irf_lin_low[, , 4]))
  figure20cseries_low <- as.data.frame(figure20cseries_low$V1)
  
  # Create Figure 20d Series 
  figure20dseries_high <- as.data.frame(t(results_lin_p70906$irf_lin_up[, , 4]))
  figure20dseries_high <- as.data.frame(figure20dseries_high$V1)
  
  figure20dseries_mean <- as.data.frame(t(results_lin_p70906$irf_lin_mean[, , 4]))
  figure20dseries_mean <- as.data.frame(figure20dseries_mean$V1)
  
  figure20dseries_low <- as.data.frame(t(results_lin_p70906$irf_lin_low[, , 4]))
  figure20dseries_low <- as.data.frame(figure20dseries_low$V1)
  
  # Create Figure 20e Series 
  figure20eseries_high <- as.data.frame(t(results_lin_p50706$irf_lin_up[, , 4]))
  figure20eseries_high <- as.data.frame(figure20eseries_high$V1)
  
  figure20eseries_mean <- as.data.frame(t(results_lin_p50706$irf_lin_mean[, , 4]))
  figure20eseries_mean <- as.data.frame(figure20eseries_mean$V1)
  
  figure20eseries_low <- as.data.frame(t(results_lin_p50706$irf_lin_low[, , 4]))
  figure20eseries_low <- as.data.frame(figure20eseries_low$V1)
  
  # Create Figure 20f Series 
  figure20fseries_high <- as.data.frame(t(results_lin_bot506$irf_lin_up[, , 4]))
  figure20fseries_high <- as.data.frame(figure20fseries_high$V1)
  
  figure20fseries_mean <- as.data.frame(t(results_lin_bot506$irf_lin_mean[, , 4]))
  figure20fseries_mean <- as.data.frame(figure20fseries_mean$V1)
  
  figure20fseries_low <- as.data.frame(t(results_lin_bot506$irf_lin_low[, , 4]))
  figure20fseries_low <- as.data.frame(figure20fseries_low$V1)
  
  # Merge all series together
  figure20series <- list(figure20aseries_high, figure20aseries_mean, figure20aseries_low, figure20bseries_high, figure20bseries_mean, figure20bseries_low, figure20cseries_high, figure20cseries_mean, figure20cseries_low, figure20dseries_high, figure20dseries_mean, figure20dseries_low, figure20eseries_high, figure20eseries_mean, figure20eseries_low, figure20fseries_high, figure20fseries_mean, figure20fseries_low)
  figure20series <- as.data.frame(figure20series)
  figure20series <- (mutate(figure20series, quarter=row_number()))
  quarter <- as.data.frame(figure20series$quarter)
  figure20series <- list(quarter, figure20aseries_high, figure20aseries_mean, figure20aseries_low, figure20bseries_high, figure20bseries_mean, figure20bseries_low, figure20cseries_high, figure20cseries_mean, figure20cseries_low, figure20dseries_high, figure20dseries_mean, figure20dseries_low, figure20eseries_high, figure20eseries_mean, figure20eseries_low, figure20fseries_high, figure20fseries_mean, figure20fseries_low)
  figure20series <- as.data.frame(figure20series)
  figure20series <- rename(figure20series, quarters = `figure20series.quarter`, figure20ahigh = `figure20aseries_high.V1`, figure20bhigh = `figure20bseries_high.V1`, figure20chigh = `figure20cseries_high.V1`, figure20dhigh = `figure20dseries_high.V1`, figure20ehigh = `figure20eseries_high.V1`, figure20fhigh = `figure20fseries_high.V1`, figure20amean = `figure20aseries_mean.V1`, figure20bmean = `figure20bseries_mean.V1`, figure20cmean = `figure20cseries_mean.V1`, figure20dmean = `figure20dseries_mean.V1`, figure20emean = `figure20eseries_mean.V1`, figure20fmean = `figure20fseries_mean.V1`, figure20alow = `figure20aseries_low.V1`, figure20blow = `figure20bseries_low.V1`, figure20clow = `figure20cseries_low.V1`, figure20dlow = `figure20dseries_low.V1`, figure20elow = `figure20eseries_low.V1`, figure20flow = `figure20fseries_low.V1`)
  
  # Export the series 
  write.csv(figure20series, "figure20series.csv")
}
export_figure20_series()

export_figure21_series <- function()
{
  # Create Figure 21a Series 
  figure21aseries_high <- as.data.frame(t(results_lin_top016$irf_lin_up[, , 5]))
  figure21aseries_high <- as.data.frame(figure21aseries_high$V1)
  
  figure21aseries_mean <- as.data.frame(t(results_lin_top016$irf_lin_mean[, , 5]))
  figure21aseries_mean <- as.data.frame(figure21aseries_mean$V1)
  
  figure21aseries_low <- as.data.frame(t(results_lin_top016$irf_lin_low[, , 5]))
  figure21aseries_low <- as.data.frame(figure21aseries_low$V1)
  
  # Create Figure 21b Series 
  figure21bseries_high <- as.data.frame(t(results_lin_top16$irf_lin_up[, , 5]))
  figure21bseries_high <- as.data.frame(figure21bseries_high$V1)
  
  figure21bseries_mean <- as.data.frame(t(results_lin_top16$irf_lin_mean[, , 5]))
  figure21bseries_mean <- as.data.frame(figure21bseries_mean$V1)
  
  figure21bseries_low <- as.data.frame(t(results_lin_top16$irf_lin_low[, , 5]))
  figure21bseries_low <- as.data.frame(figure21bseries_low$V1)
  
  # Create Figure 21c Series 
  figure21cseries_high <- as.data.frame(t(results_lin_top106$irf_lin_up[, , 5]))
  figure21cseries_high <- as.data.frame(figure21cseries_high$V1)
  
  figure21cseries_mean <- as.data.frame(t(results_lin_top106$irf_lin_mean[, , 5]))
  figure21cseries_mean <- as.data.frame(figure21cseries_mean$V1)
  
  figure21cseries_low <- as.data.frame(t(results_lin_top106$irf_lin_low[, , 5]))
  figure21cseries_low <- as.data.frame(figure21cseries_low$V1)
  
  # Create Figure 21d Series 
  figure21dseries_high <- as.data.frame(t(results_lin_p70906$irf_lin_up[, , 5]))
  figure21dseries_high <- as.data.frame(figure21dseries_high$V1)
  
  figure21dseries_mean <- as.data.frame(t(results_lin_p70906$irf_lin_mean[, , 5]))
  figure21dseries_mean <- as.data.frame(figure21dseries_mean$V1)
  
  figure21dseries_low <- as.data.frame(t(results_lin_p70906$irf_lin_low[, , 5]))
  figure21dseries_low <- as.data.frame(figure21dseries_low$V1)
  
  # Create Figure 21e Series 
  figure21eseries_high <- as.data.frame(t(results_lin_p50706$irf_lin_up[, , 5]))
  figure21eseries_high <- as.data.frame(figure21eseries_high$V1)
  
  figure21eseries_mean <- as.data.frame(t(results_lin_p50706$irf_lin_mean[, , 5]))
  figure21eseries_mean <- as.data.frame(figure21eseries_mean$V1)
  
  figure21eseries_low <- as.data.frame(t(results_lin_p50706$irf_lin_low[, , 5]))
  figure21eseries_low <- as.data.frame(figure21eseries_low$V1)
  
  # Create Figure 21f Series 
  figure21fseries_high <- as.data.frame(t(results_lin_bot506$irf_lin_up[, , 5]))
  figure21fseries_high <- as.data.frame(figure21fseries_high$V1)
  
  figure21fseries_mean <- as.data.frame(t(results_lin_bot506$irf_lin_mean[, , 5]))
  figure21fseries_mean <- as.data.frame(figure21fseries_mean$V1)
  
  figure21fseries_low <- as.data.frame(t(results_lin_bot506$irf_lin_low[, , 5]))
  figure21fseries_low <- as.data.frame(figure21fseries_low$V1)
  
  # Merge all series together
  figure21series <- list(figure21aseries_high, figure21aseries_mean, figure21aseries_low, figure21bseries_high, figure21bseries_mean, figure21bseries_low, figure21cseries_high, figure21cseries_mean, figure21cseries_low, figure21dseries_high, figure21dseries_mean, figure21dseries_low, figure21eseries_high, figure21eseries_mean, figure21eseries_low, figure21fseries_high, figure21fseries_mean, figure21fseries_low)
  figure21series <- as.data.frame(figure21series)
  figure21series <- (mutate(figure21series, quarter=row_number()))
  quarter <- as.data.frame(figure21series$quarter)
  figure21series <- list(quarter, figure21aseries_high, figure21aseries_mean, figure21aseries_low, figure21bseries_high, figure21bseries_mean, figure21bseries_low, figure21cseries_high, figure21cseries_mean, figure21cseries_low, figure21dseries_high, figure21dseries_mean, figure21dseries_low, figure21eseries_high, figure21eseries_mean, figure21eseries_low, figure21fseries_high, figure21fseries_mean, figure21fseries_low)
  figure21series <- as.data.frame(figure21series)
  figure21series <- rename(figure21series, quarters = `figure21series.quarter`, figure21ahigh = `figure21aseries_high.V1`, figure21bhigh = `figure21bseries_high.V1`, figure21chigh = `figure21cseries_high.V1`, figure21dhigh = `figure21dseries_high.V1`, figure21ehigh = `figure21eseries_high.V1`, figure21fhigh = `figure21fseries_high.V1`, figure21amean = `figure21aseries_mean.V1`, figure21bmean = `figure21bseries_mean.V1`, figure21cmean = `figure21cseries_mean.V1`, figure21dmean = `figure21dseries_mean.V1`, figure21emean = `figure21eseries_mean.V1`, figure21fmean = `figure21fseries_mean.V1`, figure21alow = `figure21aseries_low.V1`, figure21blow = `figure21bseries_low.V1`, figure21clow = `figure21cseries_low.V1`, figure21dlow = `figure21dseries_low.V1`, figure21elow = `figure21eseries_low.V1`, figure21flow = `figure21fseries_low.V1`)
  
  # Export the series 
  write.csv(figure21series, "figure21series.csv")
}
export_figure21_series()

export_figure22_series <- function()
{
  # Create Figure 22a Series 
  figure22aseries_high <- as.data.frame(t(results_lin_top016$irf_lin_up[, , 6]))
  figure22aseries_high <- as.data.frame(figure22aseries_high$V1)
  
  figure22aseries_mean <- as.data.frame(t(results_lin_top016$irf_lin_mean[, , 6]))
  figure22aseries_mean <- as.data.frame(figure22aseries_mean$V1)
  
  figure22aseries_low <- as.data.frame(t(results_lin_top016$irf_lin_low[, , 6]))
  figure22aseries_low <- as.data.frame(figure22aseries_low$V1)
  
  # Create Figure 22b Series 
  figure22bseries_high <- as.data.frame(t(results_lin_top16$irf_lin_up[, , 6]))
  figure22bseries_high <- as.data.frame(figure22bseries_high$V1)
  
  figure22bseries_mean <- as.data.frame(t(results_lin_top16$irf_lin_mean[, , 6]))
  figure22bseries_mean <- as.data.frame(figure22bseries_mean$V1)
  
  figure22bseries_low <- as.data.frame(t(results_lin_top16$irf_lin_low[, , 6]))
  figure22bseries_low <- as.data.frame(figure22bseries_low$V1)
  
  # Create Figure 22c Series 
  figure22cseries_high <- as.data.frame(t(results_lin_top106$irf_lin_up[, , 6]))
  figure22cseries_high <- as.data.frame(figure22cseries_high$V1)
  
  figure22cseries_mean <- as.data.frame(t(results_lin_top106$irf_lin_mean[, , 6]))
  figure22cseries_mean <- as.data.frame(figure22cseries_mean$V1)
  
  figure22cseries_low <- as.data.frame(t(results_lin_top106$irf_lin_low[, , 6]))
  figure22cseries_low <- as.data.frame(figure22cseries_low$V1)
  
  # Create Figure 22d Series 
  figure22dseries_high <- as.data.frame(t(results_lin_p70906$irf_lin_up[, , 6]))
  figure22dseries_high <- as.data.frame(figure22dseries_high$V1)
  
  figure22dseries_mean <- as.data.frame(t(results_lin_p70906$irf_lin_mean[, , 6]))
  figure22dseries_mean <- as.data.frame(figure22dseries_mean$V1)
  
  figure22dseries_low <- as.data.frame(t(results_lin_p70906$irf_lin_low[, , 6]))
  figure22dseries_low <- as.data.frame(figure22dseries_low$V1)
  
  # Create Figure 22e Series 
  figure22eseries_high <- as.data.frame(t(results_lin_p50706$irf_lin_up[, , 6]))
  figure22eseries_high <- as.data.frame(figure22eseries_high$V1)
  
  figure22eseries_mean <- as.data.frame(t(results_lin_p50706$irf_lin_mean[, , 6]))
  figure22eseries_mean <- as.data.frame(figure22eseries_mean$V1)
  
  figure22eseries_low <- as.data.frame(t(results_lin_p50706$irf_lin_low[, , 6]))
  figure22eseries_low <- as.data.frame(figure22eseries_low$V1)
  
  # Create Figure 22f Series 
  figure22fseries_high <- as.data.frame(t(results_lin_bot506$irf_lin_up[, , 6]))
  figure22fseries_high <- as.data.frame(figure22fseries_high$V1)
  
  figure22fseries_mean <- as.data.frame(t(results_lin_bot506$irf_lin_mean[, , 6]))
  figure22fseries_mean <- as.data.frame(figure22fseries_mean$V1)
  
  figure22fseries_low <- as.data.frame(t(results_lin_bot506$irf_lin_low[, , 6]))
  figure22fseries_low <- as.data.frame(figure22fseries_low$V1)
  
  # Merge all series together
  figure22series <- list(figure22aseries_high, figure22aseries_mean, figure22aseries_low, figure22bseries_high, figure22bseries_mean, figure22bseries_low, figure22cseries_high, figure22cseries_mean, figure22cseries_low, figure22dseries_high, figure22dseries_mean, figure22dseries_low, figure22eseries_high, figure22eseries_mean, figure22eseries_low, figure22fseries_high, figure22fseries_mean, figure22fseries_low)
  figure22series <- as.data.frame(figure22series)
  figure22series <- (mutate(figure22series, quarter=row_number()))
  quarter <- as.data.frame(figure22series$quarter)
  figure22series <- list(quarter, figure22aseries_high, figure22aseries_mean, figure22aseries_low, figure22bseries_high, figure22bseries_mean, figure22bseries_low, figure22cseries_high, figure22cseries_mean, figure22cseries_low, figure22dseries_high, figure22dseries_mean, figure22dseries_low, figure22eseries_high, figure22eseries_mean, figure22eseries_low, figure22fseries_high, figure22fseries_mean, figure22fseries_low)
  figure22series <- as.data.frame(figure22series)
  figure22series <- rename(figure22series, quarters = `figure22series.quarter`, figure22ahigh = `figure22aseries_high.V1`, figure22bhigh = `figure22bseries_high.V1`, figure22chigh = `figure22cseries_high.V1`, figure22dhigh = `figure22dseries_high.V1`, figure22ehigh = `figure22eseries_high.V1`, figure22fhigh = `figure22fseries_high.V1`, figure22amean = `figure22aseries_mean.V1`, figure22bmean = `figure22bseries_mean.V1`, figure22cmean = `figure22cseries_mean.V1`, figure22dmean = `figure22dseries_mean.V1`, figure22emean = `figure22eseries_mean.V1`, figure22fmean = `figure22fseries_mean.V1`, figure22alow = `figure22aseries_low.V1`, figure22blow = `figure22bseries_low.V1`, figure22clow = `figure22cseries_low.V1`, figure22dlow = `figure22dseries_low.V1`, figure22elow = `figure22eseries_low.V1`, figure22flow = `figure22fseries_low.V1`)
  
  # Export the series 
  write.csv(figure22series, "figure22series.csv")
}
export_figure22_series()

export_figure23_series <- function()
{
  # Create Figure 23a Series 
  figure23aseries_high <- as.data.frame(t(results_lin_top016$irf_lin_up[, , 7]))
  figure23aseries_high <- as.data.frame(figure23aseries_high$V1)
  
  figure23aseries_mean <- as.data.frame(t(results_lin_top016$irf_lin_mean[, , 7]))
  figure23aseries_mean <- as.data.frame(figure23aseries_mean$V1)
  
  figure23aseries_low <- as.data.frame(t(results_lin_top016$irf_lin_low[, , 7]))
  figure23aseries_low <- as.data.frame(figure23aseries_low$V1)
  
  # Create Figure 23b Series 
  figure23bseries_high <- as.data.frame(t(results_lin_top16$irf_lin_up[, , 7]))
  figure23bseries_high <- as.data.frame(figure23bseries_high$V1)
  
  figure23bseries_mean <- as.data.frame(t(results_lin_top16$irf_lin_mean[, , 7]))
  figure23bseries_mean <- as.data.frame(figure23bseries_mean$V1)
  
  figure23bseries_low <- as.data.frame(t(results_lin_top16$irf_lin_low[, , 7]))
  figure23bseries_low <- as.data.frame(figure23bseries_low$V1)
  
  # Create Figure 23c Series 
  figure23cseries_high <- as.data.frame(t(results_lin_top106$irf_lin_up[, , 7]))
  figure23cseries_high <- as.data.frame(figure23cseries_high$V1)
  
  figure23cseries_mean <- as.data.frame(t(results_lin_top106$irf_lin_mean[, , 7]))
  figure23cseries_mean <- as.data.frame(figure23cseries_mean$V1)
  
  figure23cseries_low <- as.data.frame(t(results_lin_top106$irf_lin_low[, , 7]))
  figure23cseries_low <- as.data.frame(figure23cseries_low$V1)
  
  # Create Figure 23d Series 
  figure23dseries_high <- as.data.frame(t(results_lin_p70906$irf_lin_up[, , 7]))
  figure23dseries_high <- as.data.frame(figure23dseries_high$V1)
  
  figure23dseries_mean <- as.data.frame(t(results_lin_p70906$irf_lin_mean[, , 7]))
  figure23dseries_mean <- as.data.frame(figure23dseries_mean$V1)
  
  figure23dseries_low <- as.data.frame(t(results_lin_p70906$irf_lin_low[, , 7]))
  figure23dseries_low <- as.data.frame(figure23dseries_low$V1)
  
  # Create Figure 23e Series 
  figure23eseries_high <- as.data.frame(t(results_lin_p50706$irf_lin_up[, , 7]))
  figure23eseries_high <- as.data.frame(figure23eseries_high$V1)
  
  figure23eseries_mean <- as.data.frame(t(results_lin_p50706$irf_lin_mean[, , 7]))
  figure23eseries_mean <- as.data.frame(figure23eseries_mean$V1)
  
  figure23eseries_low <- as.data.frame(t(results_lin_p50706$irf_lin_low[, , 7]))
  figure23eseries_low <- as.data.frame(figure23eseries_low$V1)
  
  # Create Figure 23f Series 
  figure23fseries_high <- as.data.frame(t(results_lin_bot506$irf_lin_up[, , 7]))
  figure23fseries_high <- as.data.frame(figure23fseries_high$V1)
  
  figure23fseries_mean <- as.data.frame(t(results_lin_bot506$irf_lin_mean[, , 7]))
  figure23fseries_mean <- as.data.frame(figure23fseries_mean$V1)
  
  figure23fseries_low <- as.data.frame(t(results_lin_bot506$irf_lin_low[, , 7]))
  figure23fseries_low <- as.data.frame(figure23fseries_low$V1)
  
  # Merge all series together
  figure23series <- list(figure23aseries_high, figure23aseries_mean, figure23aseries_low, figure23bseries_high, figure23bseries_mean, figure23bseries_low, figure23cseries_high, figure23cseries_mean, figure23cseries_low, figure23dseries_high, figure23dseries_mean, figure23dseries_low, figure23eseries_high, figure23eseries_mean, figure23eseries_low, figure23fseries_high, figure23fseries_mean, figure23fseries_low)
  figure23series <- as.data.frame(figure23series)
  figure23series <- (mutate(figure23series, quarter=row_number()))
  quarter <- as.data.frame(figure23series$quarter)
  figure23series <- list(quarter, figure23aseries_high, figure23aseries_mean, figure23aseries_low, figure23bseries_high, figure23bseries_mean, figure23bseries_low, figure23cseries_high, figure23cseries_mean, figure23cseries_low, figure23dseries_high, figure23dseries_mean, figure23dseries_low, figure23eseries_high, figure23eseries_mean, figure23eseries_low, figure23fseries_high, figure23fseries_mean, figure23fseries_low)
  figure23series <- as.data.frame(figure23series)
  figure23series <- rename(figure23series, quarters = `figure23series.quarter`, figure23ahigh = `figure23aseries_high.V1`, figure23bhigh = `figure23bseries_high.V1`, figure23chigh = `figure23cseries_high.V1`, figure23dhigh = `figure23dseries_high.V1`, figure23ehigh = `figure23eseries_high.V1`, figure23fhigh = `figure23fseries_high.V1`, figure23amean = `figure23aseries_mean.V1`, figure23bmean = `figure23bseries_mean.V1`, figure23cmean = `figure23cseries_mean.V1`, figure23dmean = `figure23dseries_mean.V1`, figure23emean = `figure23eseries_mean.V1`, figure23fmean = `figure23fseries_mean.V1`, figure23alow = `figure23aseries_low.V1`, figure23blow = `figure23bseries_low.V1`, figure23clow = `figure23cseries_low.V1`, figure23dlow = `figure23dseries_low.V1`, figure23elow = `figure23eseries_low.V1`, figure23flow = `figure23fseries_low.V1`)
  
  # Export the series 
  write.csv(figure23series, "figure23series.csv")
}
export_figure23_series()

#----------------------------Figure 18-23 Chart Code------------------------####

figure18 <- read_csv("[your home directory]/IRF_Figures/OutputData/figure18series.csv") 
figure19 <- read_csv("[your home directory]/IRF_Figures/OutputData/figure19series.csv") 
figure20 <- read_csv("[your home directory]/IRF_Figures/OutputData/figure20series.csv") 
figure21 <- read_csv("[your home directory]/IRF_Figures/OutputData/figure21series.csv") 
figure22 <- read_csv("[your home directory]/IRF_Figures/OutputData/figure22series.csv") 
figure23 <- read_csv("[your home directory]/IRF_Figures/OutputData/figure23series.csv") 
figure18data <- figure18 %>% filter(quarters<17)
figure19data <- figure19 %>% filter(quarters<17)
figure20data <- figure20 %>% filter(quarters<17)
figure21data <- figure21 %>% filter(quarters<17)
figure22data <- figure22 %>% filter(quarters<17)
figure23data <- figure23 %>% filter(quarters<17)
theme_if_policystyle <- function() {
  small_margins <- ggplot2::theme(plot.margin = unit(c(2, 2, 2, 3.25), "mm"))
  
  theme_if_policystyle <- ggplot2::theme_classic() +
    theme(
      #text = element_text(family = "Arial"), commented out by m1aem03 due to user reported issues 10/17/2023
      axis.ticks.length = unit(-0.2, "cm"),
      # Remove panel border
      panel.border = ggplot2::element_blank(),
      plot.background = ggplot2::element_rect(fill = "transparent", colour = NA),
      plot.margin = ggplot2::unit(c(0.04, 0.175, 0, 0.175), "in"),
      
      axis.text = ggplot2::element_text(size = 6),
      axis.text.x = ggplot2::element_text(margin = ggplot2::margin(t = 3), color = "black", size = 6), # 10
      axis.text.y = ggplot2::element_text(margin = ggplot2::margin(r = 3), color = "black", size = 6), # 9.5
      axis.text.x.bottom = ggplot2::element_text(size = 6),
      axis.text.y.left = ggplot2::element_text(margin = ggplot2::margin(r = 3), size = 6, hjust = 1), # update 1.3.1 to handle y margin spacing issues
      axis.text.y.right = ggplot2::element_text(margin = ggplot2::margin(l = 3), size = 6, hjust = 1), # 9.5
      axis.title = ggplot2::element_text(size = 6, face = "plain"),
      
      plot.title = ggplot2::element_text(hjust = .5, size = 8, face = "bold"),
      plot.caption = ggplot2::element_text(size = 6, hjust = 0),
      plot.subtitle = ggplot2::element_text(size = 6, hjust = 1, margin = margin(0, 0, 0, 0, "pt")),
      
      axis.ticks = ggplot2::element_line(linewidth = 0.3, color = "black"),
      axis.line = ggplot2::element_line(linewidth = 0.3, color = "black"),
      axis.ticks.x = ggplot2::element_line(linewidth = 0.3, color = "black"),
      axis.line.x = ggplot2::element_line(linewidth = 0.3, color = "black"),
      axis.line.y = ggplot2::element_line(linewidth = 0.3, color = "black"),
      axis.ticks.y = ggplot2::element_line(linewidth = 0.3, color = "black"),
      
      legend.text = ggplot2::element_text(size = 8),
      legend.key.size = ggplot2::unit(0.4, 'cm'),
      legend.spacing.x = unit(0.05, 'cm'),
      legend.spacing.y = unit(-0.15, 'cm'),
      legend.key = ggplot2::element_blank(), # removes the legend title
      legend.title = ggplot2::element_blank(),
      legend.background = ggplot2::element_blank()
    )
  
  return(theme_if_policystyle)
}
light_color_set <- c("#afc98d", "#729cc3", "#cec68f", "darkgrey", "orange") 

runcharts <- function(){
#------------------------------Figure18a------------------------------------####
figure18a <- ggplot(data = figure18data, aes(x = quarters)) +
    geom_ribbon(aes(y = figure18amean, ymin = figure18ahigh, ymax = figure18alow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
    geom_line(aes(y = figure18amean, color = "Impulse response"), size = .8) +
    geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
    scale_color_manual(name = "", values = c("Impulse response" = "black")) +
    scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
    ylab(NULL) +
    xlab("Quarters") +
    labs(title = "(A) 0.1% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
    scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0022, .002)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
    scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
    theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
    theme(legend.position = c(.25, 0.85))

figure18a
#------------------------------Figure18b------------------------------------####
figure18b <- ggplot(data = figure18data, aes(x = quarters)) +
    geom_ribbon(aes(y = figure18bmean, ymin = figure18bhigh, ymax = figure18blow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
    geom_line(aes(y = figure18bmean, color = "Impulse response"), size = .8) +
    geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
    scale_color_manual(name = "", values = c("Impulse response" = "black")) +
    scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
    ylab(NULL) +
    xlab("Quarters") +
    labs(title = "(B) 99-99.9% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
    scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0022, .002)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
    scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
    theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
    theme(legend.position = c(.25, 0.85))

figure18b
#------------------------------Figure18c------------------------------------####
figure18c <- ggplot(data = figure18data, aes(x = quarters)) +
    geom_ribbon(aes(y = figure18cmean, ymin = figure18chigh, ymax = figure18clow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
    geom_line(aes(y = figure18cmean, color = "Impulse response"), size = .8) +
    geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
    scale_color_manual(name = "", values = c("Impulse response" = "black")) +
    scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
    ylab(NULL) +
    xlab("Quarters") +
    labs(title = "(C) 90-99% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
    scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0022, .002)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
    scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
    theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
    theme(legend.position = c(.25, 0.85))

figure18c
#------------------------------Figure18d------------------------------------####
figure18d <- ggplot(data = figure18data, aes(x = quarters)) +
    geom_ribbon(aes(y = figure18dmean, ymin = figure18dhigh, ymax = figure18dlow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
    geom_line(aes(y = figure18dmean, color = "Impulse response"), size = .8) +
    geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
    scale_color_manual(name = "", values = c("Impulse response" = "black")) +
    scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
    ylab(NULL) +
    xlab("Quarters") +
    labs(title = "(D) 70-90% Wealth Share", subtitle = "") +
    scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0022, .002)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
    scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
    theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
    theme(legend.position = c(.25, 0.85))

figure18d
#------------------------------Figure18e------------------------------------####
figure18e <- ggplot(data = figure18data, aes(x = quarters)) +
    geom_ribbon(aes(y = figure18emean, ymin = figure18ehigh, ymax = figure18elow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
    geom_line(aes(y = figure18emean, color = "Impulse response"), size = .8) +
    geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
    scale_color_manual(name = "", values = c("Impulse response" = "black")) +
    scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
    ylab(NULL) +
    xlab("Quarters") +
    labs(title = "(E) 50-70% Wealth Share", subtitle = "") +
    scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0022, .002)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
    scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
    theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
    theme(legend.position = c(.25, 0.85))

figure18e
#------------------------------Figure18f------------------------------------####
figure18f <- ggplot(data = figure18data, aes(x = quarters)) +
    geom_ribbon(aes(y = figure18fmean, ymin = figure18fhigh, ymax = figure18flow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
    geom_line(aes(y = figure18fmean, color = "Impulse response"), size = .8) +
    geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
    scale_color_manual(name = "", values = c("Impulse response" = "black")) +
    scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
    ylab(NULL) +
    xlab("Quarters") +
    labs(title = "(F) Bottom 50% Wealth Share", subtitle = "") +
    scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0022, .002)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
    scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
    theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
    theme(legend.position = c(.3, 0.85))

figure18f

#------------------------------Figure19a------------------------------------####
figure19a <- ggplot(data = figure19data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure19amean, ymin = figure19ahigh, ymax = figure19alow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure19amean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(A) 0.1% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.004, .004, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0042, .004)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure19a
#------------------------------Figure19b------------------------------------####
figure19b <- ggplot(data = figure19data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure19bmean, ymin = figure19bhigh, ymax = figure19blow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure19bmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(B) 99-99.9% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.004, .004, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0042, .004)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure19b
#------------------------------Figure19c------------------------------------####
figure19c <- ggplot(data = figure19data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure19cmean, ymin = figure19chigh, ymax = figure19clow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure19cmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(C) 90-99% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.004, .004, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0042, .004)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure19c
#------------------------------Figure19d------------------------------------####
figure19d <- ggplot(data = figure19data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure19dmean, ymin = figure19dhigh, ymax = figure19dlow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure19dmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(D) 70-90% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.004, .004, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0042, .004)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure19d
#------------------------------Figure19e------------------------------------####
figure19e <- ggplot(data = figure19data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure19emean, ymin = figure19ehigh, ymax = figure19elow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure19emean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(E) 50-70% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.004, .004, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0042, .004)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure19e
#------------------------------Figure19f------------------------------------####
figure19f <- ggplot(data = figure19data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure19fmean, ymin = figure19fhigh, ymax = figure19flow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure19fmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(F) Bottom 50% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.004, .004, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0042, .004)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.3, 0.85))

figure19f

#------------------------------Figure20a------------------------------------####
figure20a <- ggplot(data = figure20data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure20amean, ymin = figure20ahigh, ymax = figure20alow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure20amean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(A) 0.1% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.0015, .0035, .0010), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.002, .0035)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure20a
#------------------------------Figure20b------------------------------------####
figure20b <- ggplot(data = figure20data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure20bmean, ymin = figure20bhigh, ymax = figure20blow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure20bmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(B) 99-99.9% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.0015, .0035, .0010), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.002, .0035)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure20b
#------------------------------Figure20c------------------------------------####
figure20c <- ggplot(data = figure20data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure20cmean, ymin = figure20chigh, ymax = figure20clow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure20cmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(C) 90-99% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.0015, .0035, .0010), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.002, .0035)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure20c
#------------------------------Figure20d------------------------------------####
figure20d <- ggplot(data = figure20data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure20dmean, ymin = figure20dhigh, ymax = figure20dlow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure20dmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(D) 70-90% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.0015, .0035, .0010), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.002, .0035)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure20d
#------------------------------Figure20e------------------------------------####
figure20e <- ggplot(data = figure20data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure20emean, ymin = figure20ehigh, ymax = figure20elow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure20emean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(E) 50-70% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.0015, .0035, .0010), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.002, .0035)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure20e
#------------------------------Figure20f------------------------------------####
figure20f <- ggplot(data = figure20data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure20fmean, ymin = figure20fhigh, ymax = figure20flow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure20fmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(F) Bottom 50% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.0015, .0035, .0010), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.002, .0035)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.3, 0.85))

figure20f

#------------------------------Figure21a------------------------------------####
figure21a <- ggplot(data = figure21data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure21amean, ymin = figure21ahigh, ymax = figure21alow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure21amean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(A) 0.1% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0017, .0015)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure21a
#------------------------------Figure21b------------------------------------####
figure21b <- ggplot(data = figure21data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure21bmean, ymin = figure21bhigh, ymax = figure21blow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure21bmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(B) 99-99.9% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0017, .0015)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure21b
#------------------------------Figure21c------------------------------------####
figure21c <- ggplot(data = figure21data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure21cmean, ymin = figure21chigh, ymax = figure21clow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure21cmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(C) 90-99% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0017, .0015)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure21c
#------------------------------Figure21d------------------------------------####
figure21d <- ggplot(data = figure21data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure21dmean, ymin = figure21dhigh, ymax = figure21dlow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure21dmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(D) 70-90% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0017, .0015)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure21d
#------------------------------Figure21e------------------------------------####
figure21e <- ggplot(data = figure21data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure21emean, ymin = figure21ehigh, ymax = figure21elow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure21emean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(E) 50-70% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0017, .0015)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure21e
#------------------------------Figure21f------------------------------------####
figure21f <- ggplot(data = figure21data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure21fmean, ymin = figure21fhigh, ymax = figure21flow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure21fmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(F) Bottom 50% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.002, .002, .0005), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0017, .0015)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.3, 0.85))

figure21f

#------------------------------Figure22a------------------------------------####
figure22a <- ggplot(data = figure22data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure22amean, ymin = figure22ahigh, ymax = figure22alow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure22amean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(A) 0.1% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure22a
#------------------------------Figure22b------------------------------------####
figure22b <- ggplot(data = figure22data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure22bmean, ymin = figure22bhigh, ymax = figure22blow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure22bmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(B) 99-99.9% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure22b
#------------------------------Figure22c------------------------------------####
figure22c <- ggplot(data = figure22data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure22cmean, ymin = figure22chigh, ymax = figure22clow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure22cmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(C) 90-99% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure22c
#------------------------------Figure22d------------------------------------####
figure22d <- ggplot(data = figure22data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure22dmean, ymin = figure22dhigh, ymax = figure22dlow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure22dmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(D) 70-90% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure22d
#------------------------------Figure22e------------------------------------####
figure22e <- ggplot(data = figure22data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure22emean, ymin = figure22ehigh, ymax = figure22elow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure22emean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(E) 50-70% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure22e
#------------------------------Figure22f------------------------------------####
figure22f <- ggplot(data = figure22data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure22fmean, ymin = figure22fhigh, ymax = figure22flow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure22fmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(F) Bottom 50% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.3, 0.85))

figure22f

#------------------------------Figure23a------------------------------------####
figure23a <- ggplot(data = figure23data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure23amean, ymin = figure23ahigh, ymax = figure23alow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure23amean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(A) 0.1% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure23a
#------------------------------Figure23b------------------------------------####
figure23b <- ggplot(data = figure23data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure23bmean, ymin = figure23bhigh, ymax = figure23blow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure23bmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(B) 99-99.9% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure23b
#------------------------------Figure23c------------------------------------####
figure23c <- ggplot(data = figure23data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure23cmean, ymin = figure23chigh, ymax = figure23clow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure23cmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(C) 90-99% Wealth Share", subtitle = "", caption = paste(" ", "", sep = "\n")) +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure23c
#------------------------------Figure23d------------------------------------####
figure23d <- ggplot(data = figure23data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure23dmean, ymin = figure23dhigh, ymax = figure23dlow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure23dmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(D) 70-90% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure23d
#------------------------------Figure23e------------------------------------####
figure23e <- ggplot(data = figure23data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure23emean, ymin = figure23ehigh, ymax = figure23elow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure23emean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(E) 50-70% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.25, 0.85))

figure23e
#------------------------------Figure23f------------------------------------####
figure23f <- ggplot(data = figure23data, aes(x = quarters)) +
  geom_ribbon(aes(y = figure23fmean, ymin = figure23fhigh, ymax = figure23flow, fill = "90% confidence band"), colour = "darkgrey", alpha = 0.4, size = .5, outline.type = "both") + # this creates the grey ribbon around the line
  geom_line(aes(y = figure23fmean, color = "Impulse response"), size = .8) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .3) +
  scale_color_manual(name = "", values = c("Impulse response" = "black")) +
  scale_fill_manual(name = "", values = c("90% confidence band" = light_color_set[4])) +
  ylab(NULL) +
  xlab("Quarters") +
  labs(title = "(F) Bottom 50% Wealth Share", subtitle = "") +
  scale_y_continuous(position = "right", breaks = seq(-.002, .003, .001), sec.axis = dup_axis(labels = NULL), expand = c(0, 0), limits = c(-.0025, .003)) + # Sets y axis limits and ends the y-axis with a tick mark (caps the Y-axis) 
  scale_x_continuous(expand = c(.035, .035), breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16), labels = c("","2","","4","","6","","8","","10","","12","","14","","16"))  + # Sets x axis limits
  theme_if_policystyle() + # --Theme from /if/appl/R/Functions/releases/1.4.0/theme_if_policystyle.r
  theme(legend.position = c(.3, 0.85))

figure23f

# Create Figures ####
fig18 <- ggarrange(figure18a, figure18b, figure18c, figure18d, figure18e, figure18f, ncol = 3, nrow = 2, legend = "bottom", common.legend = TRUE)
fig19 <- ggarrange(figure19a, figure19b, figure19c, figure19d, figure19e, figure19f, ncol = 3, nrow = 2, legend = "bottom", common.legend = TRUE) 
fig20 <- ggarrange(figure20a, figure20b, figure20c, figure20d, figure20e, figure20f, ncol = 3, nrow = 2, legend = "bottom", common.legend = TRUE)
fig21 <- ggarrange(figure21a, figure21b, figure21c, figure21d, figure21e, figure21f, ncol = 3, nrow = 2, legend = "bottom", common.legend = TRUE)
fig22 <- ggarrange(figure22a, figure22b, figure22c, figure22d, figure22e, figure22f, ncol = 3, nrow = 2, legend = "bottom", common.legend = TRUE)
fig23 <- ggarrange(figure23a, figure23b, figure23c, figure23d, figure23e, figure23f, ncol = 3, nrow = 2, legend = "bottom", common.legend = TRUE)

# Save Figures ####
setwd("[your home directory]/IRF_Figures/AppendixFigures/")
ggsave(fig18, file = "figure18.eps", height = 6, width = 9, units = "in", device = cairo_ps)
ggsave(fig19, file = "figure19.eps", height = 6, width = 9, units = "in", device = cairo_ps)
ggsave(fig20, file = "figure20.eps", height = 6, width = 9, units = "in", device = cairo_ps)
ggsave(fig21, file = "figure21.eps", height = 6, width = 9, units = "in", device = cairo_ps)
ggsave(fig22, file = "figure22.eps", height = 6, width = 9, units = "in", device = cairo_ps)
ggsave(fig23, file = "figure23.eps", height = 6, width = 9, units = "in", device = cairo_ps)
}
runcharts()